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Abstract 

One  of  the  major  problems  in  dealing  with  interacting  finite  populations  of  agents,  such  as 
molecules  in  chemistry  or  people  in  populations,  is  that  there  always  exists  a  probability 
of  one  species  or  state  going  extinct.  Predicting  the  probability  of  extinction  requires  a 
knowledge  of  how  the  dynamics  progresses  towards  extinction. The  path  that  optimizes  the 
probability  to  extinction  is  defined  to  be  the  optimal  path.  Here  we  present  an  algorithm  for 
constructing,  or  growing,  the  optimal  path  to  extinction  in  systems  of  arbitrary  dimensions. 
The  algorithm  relies  on  the  calculation  of  finite-time  Lyapunov  exponents  (FTLE),  which 
provide  a  quantitative  measure  of  how  sensitively  a  system’s  behavior  depends  on  initial 
conditions  [1].  We  also  present  an  efficient  method  for  approximating  the  FTLE  field  of  a 
dynamical  system. 
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1  Introduction 


Undoubtedly,  the  key  to  understanding  effective  strategies  for  expediting  the  extinction 
of  an  epidemic  lies  in  developing  an  accurate  model  of  disease  evolution  within  a  popula¬ 
tion.  Traditionally,  the  spread  of  an  epidemic  within  a  population  has  been  modeled  as  a 
stochastic  system,  where  noise  arises  from  discrete  interactions  between  individuals. 

We  will  instead  employ  an  eikonal  approximation  to  recast  the  problem  in  terms  of  an 
effective  classical  Hamiltonian  system  [1].  It  turns  out  that  there  exists  a  specific  sequence 
of  noise  that  transitions  the  system  from  an  endemic  to  extinct  state,  representing  the  most 
probable,  rare  sequence  of  such  events.  This  heteroclinic  trajectory  in  the  phase  space  of 
the  Hamiltonian  flow  is  known  as  the  optimal  path  to  extinction  (OPE). 

Analysis  of  optimal  paths  to  extinction  has  the  potential  to  reveal  new  strategies  for 
controlling  the  spread  of  disease  and  monitoring  the  progress  of  existing  protocols  [1]. 

This  report  presents  an  algorithm  for  growing  the  optimal  path  to  extinction  in  systems 
of  arbitrary  dimensions.  The  algorithm  relies  on  the  calculation  of  finite-time  Lyapunov 
exponents  (FTLE),  which  provide  a  quantitative  measure  of  how  sensitively  a  system’s 
behavior  depends  on  initial  conditions  [1].  We  will  also  present  an  efficient  method  for 
approximating  the  FTLE  field  of  a  dynamical  system. 

The  algorithm  presented  has  successfully  grown  optimal  paths  to  extinction  for  a  num¬ 
ber  of  two  dimensional  systems,  and  shows  promise  for  working  with  four  dimensional 
systems  as  well.  However,  additional  experiments  are  required  before  accuracy  in  higher 
dimensional  systems  can  be  confirmed. 

2  Finite-Time  Lyapunov  Exponent  Approximation 

Continuous  dynamical  systems  have  quantities,  known  as  Lyapunov  exponents,  which  pro¬ 
vide  a  quantitative  measure  of  how  sensitively  a  system’s  behavior  depends  on  small  per¬ 
turbations  to  initial  conditions.  Traditionally,  Lyapunov  exponents  measure  sensitivity  in 
the  infinite  time  limit;  finite  time  Lyapunov  exponents  are  simply  Lyapunov  exponents 
evaluated  on  a  finite  time  interval  [4] . 

A  detailed  exploration  of  FTLE  calculations  can  be  found  in  [1]  and  [4].  Below,  we 
present  a  method  for  approximating  FTLE  quickly  and  efficiently,  referred  to  as  the  “Poor 
Man’s  FTLE”. 

2.1  Poor  Man’s  FTLE  Algorithm 

Assume  x  and  p  are  state  variables,  representing  a  point  in  a  Hamiltonian  system’s  phase 
space. 

1.  Evaluate  equations  of  motion  starting  at  (x,p)  for  time  T  with  integration  time  step 
At,  and  record  the  final  position  as  (x',p'). 
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2.  Pick  a  random  point  at  distance  e  away  from  (x,p),  called  (xe,p€),  where  e  is  very 
small  relative  to  the  dynamics  of  the  system. 

3.  Evaluate  equations  of  motion  starting  at  (xe,pe)  for  time  T  with  integration  time 
step  At,  and  record  the  final  position  as  (x'e,pfe). 

4.  Let  S  represent  the  Euclidean  distance  between  the  two  end  points  of  the  trajectories, 
i.e.  S  =  dist  [(V,p'),  (x',p')]. 

5.  Calculate  the  FTLE  at  (x,p)  as 


1  i  5 
a  =  — -  m  - . 

\T\  e 


A  graphical  representation  of  this  algorithm  is  shown  in  Figure  1. 


Figure  1:  Graphical  representation  of  Poor  Man’s  FTLE  algorithm.  The  red  and  green 
dotted  lines  represent  the  trajectories  of  points  (x,p)  and  (xe,pe)  through  phase  space  over 
time  period  T. 


3  Growing  Optimal  Path  to  Extinction  -  The  Algorithm 

Schwartz,  et.  al.  have  shown  that  the  local  maximums  of  FTLE  fields  correspond  to  the 
optimal  path  to  extinction,  i.e.  the  optimal  path  to  extinction  is  the  most  sensitive  to  initial 
data.  In  fact,  finding  ridges  of  local  maxima  through  phase  space  will  lead  us  towards  the 
optimal  path  to  extinction,  if  one  exists  [4]. 
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Previous  work  has  focused  on  evaluating  fields  of  FTLE  values,  and  then  picking  out 
the  optimal  paths  to  extinction  based  on  local  maxima  [1].  However,  this  process  is  a 
time  consuming  and  inefficient  method  to  detect  an  optimal  path  to  extinction,  because 
it  requires  evaluating  FTLE’s  over  regions  of  phase  space  that  are  unrelated  to  the  OPE. 
Additionally,  these  FTLE  fields  are  difficult  to  visualize  in  higher  dimensions  [5]. 

Instead,  we  present  an  algorithm  to  grow  the  optimal  path  to  extinction  in  arbitrary  n 
dimensions,  knowing  only  the  endemic  and  extinct  states. 

3.1  Review  of  Terminology  and  Notation 

For  the  purposes  of  this  explanation,  we  will  use  x  to  represent  an  n-dimensional  vector, 
which  defines  a  point  in  n  —  d  +  1  dimensional  state  space.  Equivalently,  we  can  represent 
a  point  in  d  +  1  dimensional  state  space  using  hyperspherical  coordinates,  where  a  is  a 
vector  of  d  angles,  along  with  a  radius  r.  We  will  define  functions  ^(a)  and  &p(x)  to 
represent  conversions  from  polar  to  cartesian  and  cartesian  to  polar,  respectively1.  For  a 
review  of  hyperspherical  coordinates,  see  Appendix  A. 

Let  P  be  a  list  of  equally  sized  surface  partitions  of  a  d-sphere,  defined  by  intervals 
of  hyperspherical  angles2.  They  can  be  thought  of  as  bins  for  points  on  the  surface  of  a 
d-sphere.  Each  p  G  P  defines  an  equally  sized  d- volume.  Let  v  be  the  length  of  P  (i.e.  the 
number  of  partitions). 

3.2  Growing  Algorithm 

The  OPE  growing  algorithm  combines  the  fact  that  the  path  must  be  curve  along  with  the 
knowledge  that  the  local  maximum  of  the  FTLE  falls  along  the  path  to  grow  it  from  the 
endemic  state.  The  algorithm  itself  is  presented  below. 

Initialize  variables  . . . 

1.  Set  the  iterator  n  =  0. 

2.  Choose  a  starting  point  for  the  algorithm,  xo  (the  endemic  state). 

3.  Choose  an  ending  point  for  the  algorithm,  Xf  (the  extinct  state). 

4.  Choose  search  shell  radial  limits,  R.  i?max  =  is  usually  sufficient. 

Begin  procedure  . . . 

5.  Calculate  the  direction  of  the  vector  field  at  point  u. 

1For  simplicity,  these  functions  are  assumed  to  operate  only  on  points  where  \x\  =  1. 

2 Recall  that  a  d-sphere  exists  in  d+1  dimensional  space. 
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6.  Generate  S  uniformly  distributed  points  on  a  d- hemisphere,  where  the  inward  normal 
to  the  hemisphere  is  u  . 

7.  Randomly  choose  S  radii  for  the  test  points,  normally  distributed  between  Rm[n  and 

-f^max* 

8.  Remove  any  test  points  that  fall  within  a  d-sphere  of  radius  i?max  from  xK-i. 

9.  For  each3  remaining  test  point,  6,  calculate  and  store  the  FTLE  at  b  as  a using  the 
Poor  Man’s  FTLE  approximation. 

10.  Sort  the  test  points  by  their  FTLE  values,  and  extract  the  z  test  points  with  the 
largest  FTLE. 

11.  Sort  the  top  z  points  into  v  bins  based  on  spherical  angles  from  xK. 

12.  Find  the  average  angles  o:max  associated  with  the  bin  with  most  members. 

13.  Calculate  X^+x  —  i?max  ‘  ^cK^max)  T  Wn 

14.  If  \xK+i  —xf  \  >  R  max- 

•  set  K  =  K,  +  1 

•  return  to  step  6 

A  graphical  depiction  of  this  process  is  presented  in  Fig  2,  for  the  case  n  =  2. 


3  This  loop  can  be  run  in  parallel. 
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Figure  2:  Sample  progress  of  the  optimal  path  growing  algorithm  in  two  dimensional  state 
space.  Regions  of  valid  test  points  are  shaded  in  light  green  for  each  step  of  the  algorithm. 
At  each  step,  the  hemisphere  of  valid  points  has  an  inward  normal  aligned  with  the  vector 
field  at  that  point.  The  dotted  line  represents  the  grown  path  along  the  curve,  while  the 
solid  black  line  represents  the  analytical  solution. 
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4  Results 


We  now  present  the  results  of  attempting  to  grow  an  optimal  path  to  extinction  using 
the  algorithm  described  above  on  several  examples,  along  with  analytical  solutions,  if 
applicable. 

4.1  Example  1  -  Extinction  in  a  Branching- Annihilation  Process 

A  simple  system  with  intrinsic  noise  fluctuations  is  extinction  in  the  stochastic  branching- 
annihilation  process 

A  2A  and  2 A  A  0, 

where  A  and  (i  >  0  are  constant  reaction  rates  [1].  Hamilton’s  equations  are  given  as 

q  =  q[\(l  +  2p)  -  n(l+p)q] 
p  =  p[/j,(2+p)q- \(l+p)] 

where  q  is  transformation  on  the  number  of  entities,  and  p  is  the  effective  force  of  the  noise. 

It  is  simple  to  find  the  fixed  points  of  the  system  analytically  as  (A//i,  0),  (0,0)  and 
(0,  —1).  We  can  find  an  analytical  solution  for  the  optimal  path  to  extinction,  given  as 

2A(1  +  p) 
q  M2 +p) ' 

As  a  check  of  our  methods,  the  analytical  solution  for  the  optimal  path  to  extinction 
is  plotted  against  the  field  of  approximate  FTLE  values  in  Figure  3.  We  can  see  that  the 
analytical  solution  corresponds  nicely  to  the  ridge  of  local  maximum  FTLE. 

We  can  grow  this  optimal  path  without  knowing  the  analytical  solution.  Starting  at 
the  endemic  state  xo  =  (A//i,  0),  we  reach  the  extinct  state  Xf  —  (0,  —1)  by  following  the 
algorithm  above.  The  successfully  grown  curve  is  shown  in  Figure  3. 
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Figure  3:  Graphical  representation  of  running  forward  and  backwards  (averaged)  “Poor 
Man’s  FTLE”  algorithm  on  (1).  Analytical  solution  is  plotted  in  dotted  black.  Axis  are 
thickened  for  convenience. 


4.2  Example  2  -  SIS  Epidemic  Model  with  External  Fluctuations 

A  classical  Susceptilbe-Infectious-Susceptible  (SIS)  epidemiological  model  is  represented 
by  the  following  system  of  equations: 


S  =  fi  —  /jlS  +  7 1  —  PIS 
i  =  -(ii  +  i)i  +  pis 


where  fi  represents  a  constant  birth/death  rate,  ft  represents  the  contact  rate,  and  7 
represents  the  rate  of  recovery.  If  we  normalize  so  that  S  + 1  =  1,  we  can  rewrite  this  as  a 
one-dimensional  problem: 

i  =  -  (m  +  7)I  +  /3I(1-/). 

We  can  then  include  external  fluctuations  due  to  random  migrations  to/from  other  popu¬ 
lations: 


i  =  F(i)  +  n(t) 

F(-0  =  —(fi  +  7)/  +  (31(1  —  I) 

where  rj(t)  is  uncorrelated  Gaussian  noise  with  zero  mean.  Evaluating  the  Euler-Lagrange 
equation  produces  the  desired  system 

i  =  p 

P  =  (/37(1 -/)-«/)  (/5(1  -  2 :/)-«). 


where  k  —  fi  +  7  [1]. 

The  FTLE  field,  along  with  the  successfully  grown  curve,  is  shown  in  Figure  4. 

4.3  Example  3  -  SIS  Epidemic  Model  with  Internal  Fluctuations 

We  can  also  look  at  the  version  of  the  SIS  epidemic  model  where  the  form  of  the  noise  is 
not  known  beforehand.  In  this  case,  Hamilton’s  equations  become 

/=—(/, +  7)/  e~p  +  pi(l  —  I)ep 

P  =  -(m  +  l)(e~p  ~  1)  +  /3(ep  -  1)(2 1  -  1). 

No  analytical  solution  for  the  optimal  path  to  extinction  exists  for  this  system,  but 
we  can  use  our  FTLE  growing  method.  Results  obtained  from  running  the  “Poor  Man’s 
FTLE”  on  this  system,  along  with  the  grown  path,  are  presented  in  Fig.  5. 
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Grown  OPE 


Figure  4:  Graphical  representation  of  running  forward  and  backwards  (averaged)  “Poor 
Man’s  FTLE”  algorithm  on  SIS  model  with  external  fluctuations.  The  optimal  path  to 
extinction  is  easily  distinguishable. 


Figure  5:  Graphical  representation  of  running  forward  and  backwards  (averaged)  “Poor 
Man’s  FTLE”  algorithm  on  SIS  model  with  internal  fluctuations.  The  optimal  path  to 
extinction  is  easily  distinguishable. 
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4.4  Example  4  -  SIS  Epidemic  Model  without  an  Adiabatic  Assumption 

The  preceding  examples  illustrate  the  effectiveness  of  this  method  in  growing  an  optimal 
path  to  extinction  for  systems  of  two  dimensions.  We  would  now  like  to  use  this  algorithm 
to  find  an  OPE  for  a  four  dimensional  case.  We  will  use  the  SIS  model  with  internal  noise, 
but  we  will  not  make  the  adiabatic  assumption  that  S  +  I  =  1  in  order  to  remove  a  variable. 
We  begin  with  the  Hamiltonian 

H(x,p)  =  fi(ePl  —  1)  +  fix  i(e~Pl  —  1)  +  fix  2{e~P2  —  1)  +  nx2(ePl~P2  —  1)  +  /3x  iX2(eP2~Pl  —  1). 


Hamilton’s  equations  become 


x\  —  =  fiePl  —  fix \e  P1  +  K,X2ePl  p 2  —  (3x\X2eP2  Pl 

dpi 


dH 


Pi  =  =  ~lJj  (e  P1  ~  !)  “  Px<2  ( eP 2  P1  ~  !) 


x 2  —  ~px2e  p 2  —  KX2ePl  p 2  +  /3xiX2ePl  p 2 

dp2 


dH  (  -P2  l\  ( 

K  =  =  1)_K(e 


,P1  -p  2 


l)  —  /3x i  ( ep 2  Pl  —  l) 


and  the  critical  points  of  the  system  are 

/  Xl  \ 

(  1/^0  \ 

(  i  \ 

Pl 

0  1 

0 

z  = 

x2 

^endemic  — 

1 

Sa 

o 

Is 

o 

-^extinct  — 

0 

\P2  J 

V  o  J 

^ 0  / 

where  Rq  —  /?/ (fi  +  k)  . 

A  simpler  form  of  these  equations  can  be  derived  if  we  use  a  Taylor  series  expansion 
approximation  ( ex  —  x  —  1)  to  remove  the  exponentials  from  the  Hamiltonian.  We  arrive 


at: 


H(x,p)  =  ppi  (i  +  y)  -  pxipi  (i  -  y)  -  PX2P2  (l  ~  y) 


+  KX2  (Pl  -  p2)  (  1  +  ~  ~  —  )  +  fixiX2  (p2  ~  Pl)  (  1  +  P2  Pl 
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with  differentials  of 


x\  =  p,  (1  +pi)  -  nxi  (1  -  pi)  +  kx2  (pi  —P2  +  1)  -  /3xix2  (P2  —  pi  +  1) 

Pi  =  PPi  (i  -  y )  -  Px2  (P2  ~  Pi)  ^1  +  —  2  —  j 

X2  =  (P2  -  1)  -  HX2  (pi  -  P2  +  1)  +  /3xiX2  (P2  ~  Pi  +  1) 

P2  =  PP2  (f  -  y )  -  K  (pi  -  p2)  ^1  +  Pl  -  ~  /3.X i  (p2  -  pi)  ^1  +  2  Pl  ^ 

and  a  new  endemic  critical  point 


^extinct 


The  extinct  state  is  the  same  as  in  the  full  system. 

Attempts  at  running  the  optimal  path  growing  algorithm  on  the  above  systems  did  not 
terminate,  although  the  Euclidean  distance  to  the  end  point  decreased  steadily  over  time. 
A  bug  in  the  algorithm  code  may  exist,  preventing  it  from  working  in  a  higher  dimensional 
case.  Additional  review  of  the  algorithm,  as  well  as  more  control  cases  are  recommended 
to  determine  the  validity  of  this  algorithm  in  higher  dimensions. 
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A  Hyperspherical  Coordinates 


A  point  on  a  d- sphere,  which  is  embedded  in  (d  +  l)-space,  can  be  defined  either  by  d  +  1 
Cartesian  coordinates  x  =  [x\  X2  ...  Xd+ 1]  or  a  distance  from  the  origin,  i?,  along  with 
d  polar  coordinates:  ol  =  [aq  02  ...  ctd],  where  the  longitude,  oq  G  [0,  2tt)  ,  and  all  other 
colatitudes,  aq  |^>i  G  [0,7t). 

We  will  define  functions  ^(a)  and  &p(x)  to  represent  conversions  from  polar  to 
cartesian  and  cartesian  to  polar,  respectively. 

The  polar  coordinates  ol  are  defined  such  that 


d 

x\  —  R  cos  oq  sin  otj 
3= 2 
d 

x2  —  R  sin  aj 

3= 1 

d 

Xk  —  Rcosak-i  JJ  sin  oq ,  k  G  {3, . . . ,  d  +  1}. 
j=k 


To  go  the  other  direction  (cartesian  to  polar),  use  the  equations: 


,  -1  (  X2 
a  i  =  tan  — 

X\ 


otk  —  cos 


-l 


xk+ 1 


n 


■j=k+ 1  Sm  a3 


| ,  k  G  {2, ...,  d}. 


Note  that  since  the  equation  for  is  dependent  on  otk+u  if  is  most  convenient  to  solve  for 
ol  backwards;  i.e.  otd- 15  •  •  • ,  aq.  Also,  since  aq  ranges  from  0  to  2tt,  a  correction  must 

be  made  to  account  for  the  quadrant  it  is  in.  It  can  be  corrected  by 


(b  >  0  &  a  >  0)  aq 

{b  >  0  &  a  <  0)  aq  +  2tt 
{b  <  0)  oq  +  7 r 


where  a/  is  the  corrected  angle, 
instance 


R  =  \JJx~- 


The  radius,  i?,  is  analogous  to  the  three  dimensional 
x)  =  ^x\  +  ~x\  +  . . .  +  x2d  +  X2d+V 
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B  Generating  uniformly  distributed  hyperhemispheres 


To  find  a  random  point  xr  in  the  desired  region,  create  a  vector  x  of  length  n,  where 
each  entry  is  a  random  value  between  zero  and  one.  Normalize  cc,  and  multiply  by  the 
desired  radius  [3].  Repeat  as  necessary  to  generate  more  points.  This  can  be  implemented 
in  MATLAB  easily: 


S  =  1000;  %the  number  of  points  on  the  sphere  surface  we  want 
d  =  2;  %  we  want  points  on  a  d— sphere  (in  d+1  dimensions) 

C  =  randn (d+1 ,  S ) ;  %get  random  Cartesian  points 

for  ii  =  1 : S 

tempR  =  sqrt (dot (C ( : , ii) , C ( : , ii) ) ) ;  %normalize  each  point 
C  ( : , ii)  =  C  (:, ii)  . /tempR; 

end 

To  restrict  this  to  a  hyper- hemisphere,  we  need  only  limit  0  <  x\  <  1.  If  we  wish 
to  rotate  this  hyper-hemisphere  about  a  particular  direction,  we  can  just  use  the  func¬ 
tion  rotSphere(C,  a)  defined  below.  The  argument  a  represents  the  polar  angles  of  the 
hemisphere  normal  vector,  u. 


function  [Cp  R]  =  rotSphere  (C, a) 

%ROTSPHERE  Rotate  hypersphere  with  points  "C"  about  angles  defined  in  "a" 

Cp  =  zeros  (size (C) ) ; 
d  =  length (a) ; 

R  =  eye (d+1 ) ; 

for  jj  =  l:d  %create  the  composite  rotation  matrix 
R  =  R*makeRotMat (a ( j j ) ,  j j , d) ; 

end 

for  ii  =  l:size(C,2)  %for  each  point,  rotate 
Cp ( : , ii)  =  R*C ( :  ,  ii)  ; 

end 

function  R  =  makeRotMat (ang, this_dim, t ot al_d) 

%make  a  single  rotation  matrix 
R  =  eye ( total_d+l )  ; 

Rsub  =  [cos (ang)  —sin (ang); 

sin  (ang)  cos  (ang)  ] ; 

R ( this_dim : this_dim+l , this_dim : this_dim+l )  =  Rsub; 


15 


C  Hypersphere  Partitioning  and  Binning 


The  optimal  path  growing  algorithm  is  assumes  a  method  exists  for  partitioning  the  surface 
of  a  d-sphere  and  sorted  points  on  the  surface  into  such  partitions  [2].  In  the  case  d  —  1 
this  is  trivial;  one  can  simply  divide  the  polar  coordinate  0  of  a  circle  into  linearly  spaced 
partitions.  For  instance,  doing  theta  =  linspace  (0, 2*pi,  nBins)  in  MATLAB  would 
results  in  a  list  of  boundaries  of  a  circumference  partitioned  into  nBins  sections. 

Unfortunately,  this  is  not  so  simple  in  cases  where  d  >  1.  An  attempt  to  use  this 
method  in  the  case  d  —  2  is  shown  in  Fig.  6.  The  two  polar  coordinates  9  and  (j)  are 
linearly  spaced,  however  the  surface  elements  resulting  from  these  spacing  are  clearly  of 
unequal  size. 

0.999 
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0.333 

0.111 

-0.111 

-0.333 

-0.555 

-0.777 


Figure  6:  An  example  2-sphere  (in  3-dimensions)  with  linearly  spaced  polar  angles.  Linearly 
spaced  polar  angles  lead  to  inconsistent  surface  area  sizes  when  d  >  1. 

Instead,  we  need  a  method  for  partitioning  the  surface  of  a  d-sphere  into  equally  sized 
d-volume  elements.  Fortunately,  a  MATLAB  package  has  been  developed  for  this  very 
purpose.  Once  the  package  is  installed,  simply  use  the  command 
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regions  =  eq_regions_polar (d,nBins) ;  %polar  coordinate  intervals 
r.centers  =  eq_region_centers (d,nBins) ;  %polar  coordinate  centers 

to  produce  a  matrix  which  represents  nBins  equally  sized  partitions  of  a  d-sphere.  A 
2-sphere  generated  with  this  package  is  presented  in  Fig.  7. 


Figure  7:  (a)  An  example  2-sphere  with  10  equally  sized  surface  partitions,  (b)  A  stere¬ 
ographic  projection  of  a  3-sphere  divided  into  10  equally  sized  surface  partitions.  Gener¬ 
ated  by  the  Recursive  Zonal  Equal  Area  Sphere  Partitioning  Toolbox  with  the  command 
pro  ject_s3_part  it  ion  (10)  [2]. 

Once  the  partition  polar  coordinate  intervals  are  known,  sorted  points  on  the  surface 
of  a  d-sphere  into  bins  is  a  trivial  task.  For  example,  to  bin  points  on  the  surface  of  a 
3-sphere,  simply  iterate  through  each  partition  until  a  points  polar  angles  are  within  the 
of  a  partition’s  three  intervals. 
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